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We investigate dynamics in two-component Bose-Einstein condensates in the context of coupled 
Gross-Pitaevskii equations and derive results for the evolution of the total density fluctuations. Using 
these results, we show how, in many cases of interest, the dynamics can be accurately described with 
an effective one-component Gross-Pitaevskii equation for one of the components, with the trap and 
interaction coefficients determined by the relative differences in the scattering lengths. We discuss 
the model in various regimes, where it predicts breathing excitations, and the formation of vector 
solitons. An effective nonlinear evolution is predicted for some cases of current experimental interest. 
We then apply the model to construct quasi-stationary states of two-component condensates. 
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I. INTRODUCTION 

One of the most intriguing aspects of recently pro- 
duced atomic Bose-Einstein condensates p| is the ease 
with which two (or more) internal levels can be popu- 
lated, in effect forming a spin- 1/2 (or spin- 1, etc.) sys- 
tem. In particular, phase separation HyL vector soli- 
tons Q, two-component vortic es IE EL LlB , spin waves 
0, breathe-together solutions llcfTllT a nd more general 
issues regarding the dynamics |lfl lid Hi^ and stabil- 
ity Q in these systems have recently been investigated. 
Two-component BECs hold promise for a number of ap- 
plications, including coherent storage and processing of 
optical fields using stopped light techniques [Tr^|. 

The dynamics of two-component BECs is in general a 
difficult problem to study analytically, and much of the 
recent activity has focused on particular situations and 
thus made specific assumptions. For example, the rela- 
tive interaction strengths (via the binary collision scat- 
tering lengths) for inter- and intra-component collisions 
are sometimes presumed to be equal, whereas in practice 
these scattering lengths vary and depend on the atom and 
the particular levels utilized. Furthermore, these can be 
experimentally varied with Feshbach resonances |17| . In 
other studies it is assumed there is little or no spatial 
variation in the relative density between the two compo- 
nents, restricting the applicability of these calculations 
to many cases of interest. For example, in stopped light 
experiments, a highly spatially dependent superposition 
of the two internal atomic states are generated by light 
pulses, and the dynamics of this inhomogenous superpo- 
sition is essential in determining the information storage 
and processing capabilities of the system [15j . 

The need for a simple understanding of dynamics in 
two-component BECs free of such assumptions is the 
motivation for this work. We specifically consider the 
experimentally relevant situation of a single-component 
BEC (in an internal state labelled |1)) in its ground state, 
whose constituent atoms are suddenly put into spatially 
dependent superposition of two internal levels (|1) and 



2)), with the |2) condensate occupying a region inside 
the larger |1) condensate. Such a situation occurs with 
stopped light pulses or spatially dependent Raman 
pulses. When such processes occur fast compared to the 
atomic dynamics (millisecond timescales), the resulting 
two-component BEC will initially have the same density 
profile as the original single-component BEC. However, 
this is not a stationary state and will evolve. The evolu- 
tion can be modelled with coupled Gross-Pitaevskii (GP) 
equations, where each component is described with a sin- 
gle macroscopic wave functions ipi and ip2 ■ 

We find, remarkably, that the evolution of ip2 can of- 
ten be described with an effective single-component GP 
equation, with the trapping potential and interaction 
coefficients renormalized by the fractional difference in 
scattering lengths. Those these differences are generally 
small in practice, they end up dominating the motion of 
the wave functions. Depending on their sign and mag- 
nitude, our model predicts both trapping and repulsive 
effective potentials, as well as positive (repelling) or neg- 
ative (attractive) effective scattering lengths for the |2) 
condensate. 

The reduction to a single-component picture is accom- 
plished by observing that the fluctuations of the total 
density are smaller and more quickly varying than the 
evolution of each individual component. We find an 
equation of motion for the total density fluctuations and 
find that, within certain limits, these fluctuations are de- 
scribed with a simple expression which we can plug into 
the equation of motion for ip2 and derive our effective 
one-component description. While we will perform our 
calculations with one-dimensional equations for compu- 
tational simplicity, the results are equally valid and ap- 
plicable in full three-dimensional geometries. 

Using our model, much of the vast existing literature 
on single-component BEC dynamics can be easily be ap- 
plied to predict analogous behavior in two-component 
systems. To demonstrate this applicability, and test the 
accuracy of our model, we present calculations in sev- 
eral parameter regimes. For certain relative scattering 
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lengths we get a repulsive effective nonlinearity in a trap. 
In this case a simple breathing motion occurs and the ex- 
isting analytic results on collective excitations in single- 
component BECs can be applied. In other cases, we 
get an attractive nonlinearity, in which case phase sep- 
aration occurs and vector solitons form and propagate. 
Our model again allows a mapping onto existing litera- 
ture on the formation of soliton trains in single compo- 
nent BECs 0, |2(j which should improve abilities to de- 
sign and analyze experiments to observe vector solitons, 
which have hitherto not been ovserved. We also note 
that the levels |F = 1, M F = -1) and |F = 2, Mp = +1) 
in 87 Rb, used in many present day experiments [EL IT^. 
give rise to a vanishing effective nonlinearity, allowing us 
predict the evolution with a linear Schroedinger equa- 
tion. This should prove especially powerful in two and 
three dimensional cases, where it is computationally ex- 
pensive to solve the full nonlinear differential equations. 
Finally, we use our model to predict the existence of so- 
lutions whereby two overlapping condensates (each with 
arbitrary number) remain nearly stationary for very long 
times. These solutions can be thought of as generali zing 
previously discussed breathe-together solutions 0, Hlf . 
These solutions are relevant to observing spin squeezing 
[2ll | and coherent optical storage [13] in these systems. 



II. DERIVATION OF THE ONE-COMPONENT 
MODEL 

A. Description of the system 

The coupled GP equations governing the evolution are 
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where m is the mass of the atoms and the harmonic ex- 
ternal trapping potential V(z) = \mui 2 z 2 is assumed to 
be equal for each component (which is true in particular 
sub-levels for magnetic traps [12] for and for all sub-levels 
in far off-resonant optical traps [22|1. Atom-atom inter- 
actions are characterized by the C/y = 4irNh 2 aij /mA, 
where TV is the total number of condensate atoms, ay- 
are the s-wave scattering lengths for binary collisions and 
between atoms in internal states \i) and \ We are only 
accounting for dynamics in the the z dimension, which 
is valid in an elongated trap geometry where the signif- 
icant dynamics occur primarily in this dimension [23j, 
and so choose an effective transverse area A to give the 
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FIG. 1: (a) The densities N\t/j2\ 2 (solid curve, red online) 
and | 2 (dashed curve, blue online), and total density 

Np (dotted, black curve) at the times indicated. The initial 
total density po corresponds to the ground state of a pure |1) 
condensate. The dot-dashed curve shows the fluctuations of 
the total density Sp scaled by —1/S C . (b) The time evolution 
of the density profile \ip2 1 2 in this case, (c) The time evolution 
of 8p. 



correct nonlinearity . We will show how the relative dif- 
ference in the scattering lengths are key parameters in 
the evolution and so define S c = {a\ 2 — an) /an and 
$2 = {a 22 - an) /an. 

For concreteness, we consider a condensate with N — 
2.0 x 10 6 atoms in a trap with frequency io z = (2ir) 21 Hz. 
The area A = ir(4.2 /im) 2 is chosen such that the initial 
density in the center corresponds to a ground state BEC 
in a cylindrically symmetric trap with a transverse fre- 
quency uj r = 8uj z . The density profile of this initial state, 
labelled p , is plotted as the dotted curve in the first panel 
of Fig.^a). We then assume the BEC is put into a spa- 
tial superposition of two states with a Gaussian-shaped 
and slightly off-center wavefunction ip2, inside the larger 
|1) condensate [again see the first panel of Fig.QJa)]. The 
total density still given by po but, because of the spin ex- 
citation, this is no longer a stationary state. 

The resulting two-component evolution is presented 
in Fig. nja). Here we have chosen 62 — —0.03 and 
S c = —0.03. As we observe in the figure, each compo- 
nent undergoes non-trivial evolution of its density profile 
while the overall density profile (which we define to be 
p = \ijji\ 2 + \ip2 1 2 ) undergoes only very small fluctuations. 
To get a sense of the shape and magnitude of these small 
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fluctuations, the dot-dashed curves show the deviations 
from the original density profile, Sp = p — p , blown up 
by the factor —1/S C ~ 33 (the reason for this particular 
case will become apparent later). One sees that at 15 ms 
a slight fluctuation in the total density has appeared in 
the region occupied by \2) as well as near the condensate 
edges (in such a way that the total atom number is con- 
served). Figure gives an overview of the evolution 
of the density profile \ip2\ 2 and one sees that there is a 
breathing like behavior, with the width becoming larger 
and then smaller, as well as a small dipolc sloshing due 
to its original offset from the trap center. Figure [He) 
then shows the evolution of the density fluctuations Sp 
and one observes two important features. First there is 
a pair of density perturbations travelling back and forth 
across the |1) BEC (and reflecting at its boundaries), 
which, as we will see below, are travelling at the sound 
speed determined by the total BEC density. In the pan- 
els of Fig. Qfa), showing 15 ms and 45 ms, these waves 
are near the BEC boundary whereas at 30 ms they are 
crossing in the region occupied by V^l 2 - Second, there is 
a part which appears to be closely mimicking the much 
more slowly evolving profile of \ip2\ 2 (as also seen in the 
Fig. da) plots). 
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We next linearize in the velocity field v c and density 
fluctuations Sp — p — po. In this context, po is defined as 
the single component (/ = 0) stationary solution: 
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B. Calculation of density fluctuations 



Performing the linearization and eliminating v c yields a 
second order equation for the density fluctuations: 



We now analytically investigate the coupled GP equa- 
tions illll'l) to learn how this particular pattern for the 
total density fluctuations arises. Our strategy will be to 
eliminate the wavefunction ipi in favor of a hydrodynamic 
description of the total density p and total velocity field 
v c — (vi | f/^i | 2 + v 2 \ip 2 \ 2 )/p, where Vi = (fi/m)^ are the 
velocity fields of each component (the denote the gra- 
dients of the phases of the wavefunctions ipi ; we will use 
the prime symbol to indicate d/dz). We then linearize in 
the small density fluctuations Sp and consequently small 
velocity field v c . Our observation that Sp evolves on a 
fast time scale relative to ip 2 supplies an obvious separa- 
tion of time scales in the problem and allows us to then 
solve for the evolution Sp assuming ip 2 is static on this 
time scale. 

For convenience we define the relative density in \2) to 
be / = \ip 2 \ 2 /p, so 1 — / = \ipi\ 2 /p. After some lengthy, 
but straightforward, algebra one can show that Eqs. (11121) 
imply 
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where ^MFii e MF2i e KEi an< ^ e SP anc ^ are obtained by 
making the replacement p — > p in the corresponding ex- 
pressions Ij5l8(l . In this expression we have dropped terms 
involving the product of the fluctuations Sp with kinetic 
energy terms and relative scattering length differences 
~ S c , S 2 - 

Equation I|1U|) is quite widely applicable, however it 
turns out that in practice we can simplify things further; 
in the Thomas-Fermi regime the spatial derivatives of the 
background density po are small compared to the mean 
field and as a result the second line of Eq. lfTU|l can be 
neglected relative to the first line. The result is an in- 
tuitively simple picture: The first term implies the den- 
sity fluctuations obey a phonon-like dispersion term with 
usual sound speed in the condensate cq = \/Uiipo/m, 
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while the remaining terms of the first line provide var- 
ious source terms which seed these density fluctuations 
and are non-zero only in locations where there is a spin- 
excitation (that is, where / ^ 0). These source terms 
occur both because of differences in the mean field inter- 
action between the components as well as kinetic energy 
and quantum pressure in the two wavefunctions. Because 
the spin dynamics are much slower than the sound speed 
propagation, the source terms appear to be nearly sta- 
tionary to the phonons. 

Armed with this picture, we can now interpret 
Fig. |T[V,). The initial spin excitation gives non-zero 
source terms that generate phonons, which then prop- 
agate through the BEC towards the boundary. Mean- 
while, in the region of the spin excitation, the fluctuations 
are driven into a quasi-steady state solution. In an infi- 
nite medium the phonons would continue, however here 
they reflect off the boundaries and so repeatedly cross the 
area of the spin excitation. In the Fig. ^ such a cross- 
ing occurs at 30 ms, while reflections off the boundaries 
occur at 15 ms and 45 ms. 

In many cases of interest, the mean field contribu- 
tion €mfi dominates the kinetic energy contributions in 
Eq. (|10() . In this case it is easy to solve for the quasi- 
steady state solution Sp^ by setting Sp = 0: 



Sp (ss) = -pa 



(11) 



When / <C 1 this reduces to Sp^ = — <5 C |^2| 2 , which 
predicts fluctuations which are directly proportional to 
the density in |2). This approximate solution holds fairly 
well in the example, as shown by the dot-dashed curves 
in Fig. nfa). At 15 ms and 45 ms the sound waves are 
primarily at the BEC edge and this quasi-steady state 
solution dominates in the area of the spin-excitation. At 
30 ms, the sound waves are passing through the spin ex- 
citation region and are comparable to the quasi-steady 
state part. These sound waves end up having virtually 
no impact on the evolution of ip2 due to the fact that 
their evolution is on a completely different and indepen- 
dent time scale. Stated another way, though the sound 
waves overlap -02 each time they cross the BEC, the con- 
tributions from each crossing tend to be out of phase and 
their net contribution to the evolution of ip2 washes out 
to nearly zero. On the other hand the quasi-steady state 
part, given by Eq. (|1 1|) . can have a large impact on the 
evolution of ip2- 

Explicit numerical calculation of the various source 
terms reveals that the mean field term indeed dominates 
by more than an order of magnitude, meaning Eq. 1|11|) 
should hold. Furthermore, inspection of the figure shows 
that the initial peak relative density is / ~ 1/5 and the 
first term should dominate. 



C. Effective one-component GP equation 

The existence of this simple solution allows us to then 
solve for the much slower evolution of ip2 ■ Re- writing @ 
in terms of S c , 62, eliminating Vi \ 2 in favor of p and using 
our solution of © for the initial density profile po allows 
us to write: 
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where in the last line we have substituted our result for 
the quasi-steady state Sp w —5 c \ip2\ 2 - This is just a one- 
component GP equation governing ip2 with the trap and 
interaction coefficients renormalized by the 62 , S c and is 
the central result of our paper. Note that both the linear 
and nonlinear terms can have either sign, meaning that 
the qualitative behavior of i/> 2 is quite sensitive to the 
sign of the relative scattering lengths. Note also that 
the equation is consistent with our assumption that the 
ip2 dynamics are slower than the phonon dynamics when 
|*a|,|*d«l. 

We made a number of assumptions in deriving Eq. I|13|) 
and so we now turn our attention to studying the quanti- 
tative accuracy of this equation in several cases. Simulta- 
neously, this will allow us to explore a variety of qualita- 
tive distinct regimes it predicts. Addressing first our ex- 
ample from Fig. we plot in Fig.|2Ja) the evolution of the 
density profile \ip2\ 2 as predicted by Eq. (fH5|l . One sees 
no visible deviations from the full two component calcula- 
tion in Fig. mb). Note that because 6 C = $2 = —0.03 we 
get an effective repulsive interaction coefficient 0.03Un 
and trapping potential 0.03I^(z) giving rise to the breath- 
ing and dipole motion observed in the calculation. 

A quantitative comparison can be made by calculating 
the energy functionals: 



Ek = J dzW 2 \ 2 ; 

E v = S c [ dzV(z)\^\ 



Eint = -(S 2 -25 c ) J dz\ip 2 t- 
E = Ek + Ey + Eint; 



(14) 



which we plot in Fig. Efb) for both for the effective one- 
component model (thin, green curves) and the full two- 
component prediction (thick, orange curves). One sees 
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FIG. 2: (a) Evolution of \ip2\ 2 according to the effective one- 
component GP equation (11311 . (b) Contributions to the en- 
ergies, from Eq. I14|l . Ek (dotted curves), Ey (dashed), i?i n t 
(dot-dashed) and total E (solid) for the true two-component 
evolution according to Eqs. (MI2t (thick, orange curves) and 
effective one-component model Eq. I13H (thin, green curves). 
All energies are normalized by the initial total energy Eo = 
E(t). The dotted lines indicate the maximum and minimum 
values of E reached in the true two-component evolution. 
Their difference defines AE. (c) The relative fluctuations 
of E in the two-component model, normalized to the initial 
energy Eo, with the initial peak amplitude / varied. The cir- 
cled point indicates the case plotted in (a),(b). (d) The rela- 
tive energy fluctuations AE/Eo with parameters the same as 
Fig-Q but varying — 8 C and keeping 82 — S c . 

the curves track each other very closely and, furthermore, 
that all three contributions Ek, Ey , E- lnt are playing an 
important role in the evolution, indicating the evolution 
of ip2 is nonlinear in this case. The main deviation one 
observes is a very slight difference in the time scale for the 
oscillatory motion in the two cases. The total effective 
energy E is necessarily conserved for the one-component 
model, while this quantity will only be conserved for the 
true two-component calculation when Eq. I|13|) is provid- 
ing an accurate description of the evolution. Thus the 
fluctuations of this quantity (the solid, orange curve) is 
a good measure of the validity of the model and in the 
case there once sees these fluctuations AE are at a few 
percent level relative to the initial total energy Eq. In 
general, this error generally provides an estimate of the 
error in the oscillatory time scale. 

III. APPLICATIONS OF THE MODEL 

A. Prediction of breathing motion 

According to our model, Eq. 1131 . relative scattering 
lengths in the regime in Fig. (S c < and 82 — 2<5 C > 0) 
will give rise to an evolution of ^2 analogous to a trapped 
interacting single component BEC. To test this model 
and investigate the range of its applicability, in Fig.JSJc) 



we plot the effective energy fluctuations for a number of 
cases with the same parameters as Fig.|2Ia-b) but varying 
the initial peak relative density / = \ip2\ 2 /p- One sees 
that even up to / ks 0.7 the energy fluctuations are still 
less than 10 %. In Fig.|3fd) we show a series varying the 
magnitudes of 82,8c (keeping 82 — 8 C < 0) and see an 
approximate linear dependence. 

We note that because this regime generally leads to 
smooth breathing behavior, it may be well suited to per- 
forming controlled information processing via the two- 
component dynamics. In particular, the linear evolu- 
tion can be predicted by decomposing the wavefunction 
into the harmonic oscillator states of the effective poten- 
tial, and amount of nonlinearity can be controlled by 
the number of atoms coupled into |2). Furthermore, 
much can be borrowed from the vast existing litera- 
ture on one-component BEC dynamics, including ana- 
lytic predictions for the ground states in the Thomas- 
Fermi regime |'2 lj and the spectrum of linear excitations 
from this ground state [T^ . 



B. Phase separation and vector solitons 

One particularly interesting behavior of two species 
BECs is phase-separation 0, Q which is predicted to oc- 
cur when > Ui\U22- In this case it is energetically fa- 
vorable for the two components to separate into a series of 
non-overlapping domains. The manner in which this sep- 
aration takes place dynamically has not been addressed 
in detail. According to our model Eq. (|13fl . a species |2) 
contained in a condensate of another species 1 1) will effec- 
tively act as an attractive BEC when <?2 — 2<5 C < 0, which, 
to first order in |5 2 |, \ 5 C \ is equivalent to the above phase 
separation criteria. Such a case (62 = —0.09, 82 = —0.03) 
is presented in Figs.|2Ia-b). One sees that, indeed, the |2) 
condensate acts as if it has attractive interactions which 
leads to phase separation. The wavefunction tp2 first col- 
lapses then suddenly splits into two distinct soliton-like 
structures. Fig.[3Jb) shows that this evolution continues, 
with the the number of distinct structures alternating be- 
tween 1,2 and 3. The effective one-component model here 
allows us to map this two-component problem onto the 
problem of soliton train formation in single component, 
attractive interaction condensates, studied both experi- 
mentally [I3 and theoretically |2J|, and thus acts as an 
intuitively simple and quantitatively useful guide in pre- 
dicting the formation and dynamics of two-component 
(or vector) solitons The effective one-component pre- 
diction is plotted in Fig.[3{c) and one sees the model gives 
a qualitatively accurate prediction of the behavior. 

One sees in Fig. |3Jd) that the quantitative error in 
the evolution is slightly larger than in the effective repul- 
sive case with visible differences in the magnitude of the 
initial kinetic and interaction energy oscillations and es- 
pecially in the magnitude and timing of later oscillations. 
This is primarily because the relative density / grows to 
be quite large during the phase separation. The small 
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FIG. 3: Numerical simulations with S c = —0.03, 82 = —0.09, 
giving an effective attractive interaction. All conventions are 
the same as in Figs. 11121 (al Snapshots of the two component's 
densities profiles and the density fluctuations, (b) Time evo- 
lution of the density profile |"!/>2| 2 , (c) Effective one-component 
prediction (d) Energies, (e) Total energy fluctuations (nor- 
malized by the maximum value reached by E^). The circled 
point again indicates the case shown, and the other points 
indicate simulations with the same parameters, but varying 
the initial peak amplitude of /. 



/ prediction Sp = — S c \ip2\ 2 still provides a remarkably 
good estimate, as shown in Fig. Efa). While Eq. ltTT|) 
would seem to imply that including a next order nonlin- 
earity term ~ f 2 in the evolution could further improve 
the model, we found numerically that this did not signifi- 
cantly improve the quantitative accuracy. A further com- 
plication is that (the linearized versions of) the kinetic 
energy and quantum pressure source terms Eqs. (??-|SJ 
can be more important (relative to the mean field term 
Eqs. (??)). While explicit calculation in the case shown 
reveals they were still smaller than the mean field term by 
a factor ~ 10, this is enough to have some effect. Fig.[3lc) 
shows the energy fluctuation error (this time relative to 
the peak kinetic energy ^ max - ) since the initial total en- 
ergy Eq, being the sum of two large numbers of opposite 
sign, is near zero). One sees an approximately linear de- 



pendence of the energy fluctuations with the initial peak 
/, which saturates when the initial peak value reaches 
about 0.25. For simulations with / higher than the case 
plotted in Figs. EJa-c), the number of solitons predicted 
was incorrect. 



C. Effective repulsive potentials 

In the examples we have shown so far, we have chosen 
S c < 0, since this leads to an effective trapping potential 
for ip2- The model is accurate for the opposite case as 
well, however it predicts (and we indeed observe) that 
if>2 is pushed to the edge of the BEC due to the effec- 
tive repulsive potential. This occurs on a time scale of 
about 20 ms for 5 C = 0.03 and the initial BEC param- 
eters used in this paper. The model eventually breaks 
down when -02 reaches the BEC boundary, since the ki- 
netic energy becomes comparable to the mean field po- 
tential Unpo, and the Thomas- Fermi assumption is no 
longer valid. Similar approaches could be constructed to 
account for this boundary region in certain cases (for ex- 
ample, in 0| this was done for a case with a negligible 
nonlinearity term) but this is beyond our scope here. 



D. Vanishing nonlinearity in 87 Rb 

A particular case of interest in present day experi- 
ments is the hyperfine levels |1) = \F = l,Mp = — 1) 
and |2) = \F = 2,M F = +1} of 87 Rb. These two lev- 
els are approximately equally trapped ma gnet ically and 
have a very small inelastic collision rate [23, allowing 
them to overlap for very long times and maintain their 
coherence. Interestingly, the scattering lengths in this 
case (5 C = —0.03, 5 2 = —0.06) j2^ are such that the 
effective nonlinearity vanishes and we get a very simple 
linear evolution of "02 which could be predicted by simply 
decomposing the initial state into the eigenstates of the 
effective harmonic oscillator potential — S c V(z). There 
will be higher order terms ~ |^2| 2 , | <5 C | 2 which eventually 
introduce nonlinearity at longer times (and in fact this 
system is eventually phase separating [12| ). However the 
linear model could prove to be a powerful tool for predict- 
ing otherwise computationally expensive two-component 
evolution in two and three dimensions. 



E. Quasi-stationary solutions 

As a final application of our model, we consider the sta- 
tionary states of the effective one-component model. In a 
case with only a small density \ip2\ 2 the eigenstates of the 
harmonic oscillator potential will be stationary states of 
tp 2 . The ground state is simply a Gaussian wavefunction 
which could easily be created with a stopped light pulse 
|15| . Our model would then predict the state written 
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FIG. 4: (a) The ground state density profile |"!/>2| 2 obtained 
by finding the effective ground state of Eq. 1)13^ for a case 
with N 2 = 0.38iV, with \ipi\ 2 chosen such that the total 
density p corresponds to the original ground state profile 
po- (b) The deviations from the initial density \ip 2 {z, t) 2 — 
1^2(2, 0)| 2 are quasi-periodic, (c) The deviation parameter 
A 2 = / dz(\i> 2 (z,t)\ 2 - |0(z,O)| 2 )/ J dz\i/j 2 (z,t)\ 4 as versus 
time, (d) The maximum value reached by A2 as a function 
of the fraction of atoms in |2), N2/N. 



onto the wavefunction ip 2 is then stationary and there- 
fore quite robust with respect to storage of information. 
We have indeed observed numerically that the Gaussian 
ground state of the effective potential is stationary. Such 
an approach was used in to predict robust storage of 
optical vortex states in BECs in three-dimensional ge- 
ometries. In this work, the stability was checked with a 
full calculation of the two-component evolution, but the 
one-component model provided a way to construct an 
accurate estimate for a stable configuration. 

One of the motivations for stationary two-component 
states is the possibility of atom-atom interaction induced 
spin-squeezing, as proposed in plj . In that proposal, for 
23 Na, 5 C = —0.08, S 2 = 0. To achieve spin-squeezing, 
one must have two components, each with a compara- 
ble density, overlap for a considerable time. Ideally, the 
mean-field dynamics should be kept to a minimum to 
prevent them from washing out the squeezing dynamics. 
This latter requirement could be accomplished by choos- 
ing the number of atoms in each component consistent 
with the breathe-together solutions, whereby the rela- 
tive density between the two components (/) is constant 
across the BEC 0, 0] , while the overall density profile 
breathes. However, our model predicts quasi-stationary 
solutions for any number of atoms in |2). Furthermore, if 
one wished to write the squeezed state via slow light tech- 
niques, an inhomogenous configuration, where tp 2 is con- 
tained in tpi and vanishes at the BEC boundaries, would 
be essential to prevent spontaneous emission events. To 
our knowledge, two species BEC stationary states, with 



an inhomogenous relative density profile /, have not been 
previously predicted or experimentally observed. 

As an example, in Fig. 01a) we show a ground state of 
the effective one-component CP equation (|I3[) . obtained 
by propagating the equation in imaginary time, and hold- 
ing the number of atoms in |2) N 2 — N J dz\?p 2 \ 2 con- 
stant. In this case, the nonlinear term E- m t is quite sig- 
nificant (it dominates the kinetic energy Ek in the ground 
state). The density l^il 2 is then chosen so that the to- 
tal density po corresponds the ground state of a pure |1) 
condensate po. 

A time-dependent evolution of the full two-component 
Eqs. (11I2II with this initial state then reveals that indeed 
this configuration is nearly stationary. The only observed 
motion is a small in-phase breathing of both components, 
with a magnitude governed by 5 C . Figure Bib) shows the 
variations of the density 1021 2 from its initial value as a 
function of time, which exhibits this breathing motion. 
Thus, not only is the motion of ^2 small, but the small de- 
viations are approximately periodic and so |V>2| 2 remains 
near its initial value for very long times. The parame- 
ter A 2 , plotted in Fig. ^c) and defined in the caption, 
characterizes the total deviation and is seen to be nearly 
periodic with a maximum amplitude of about 0.025. 

We performed similar simulations for a variety of val- 
ues N2 and plotted the maximum A2 reached in each 
case. The results are presented in Fig.QJd). In the limit 
of small N 2 , where the nonlincarity is negligible the effec- 
tive ground state is truly stationary. As the nonlinearity 
becomes more important, the ground state grows due to 
effective repulsive interactions, the small breathing mo- 
tion occurs and we get the deviations A 2 . The magnitude 
of A2 saturates at around 0.025 (see Fig.^Jd)). The sat- 
uration occurs at the point the Thomas-Fermi solution 
[2^ of the effective one-component model becomes accu- 
rate. The saturation value should scale roughly linearly 
with \8 C \. 

In the regime N 2 ^ Ni, the effective ground state even- 
tually becomes larger than the original condensate, in 
which case the present approach fails. However, in that 
limit, the solutions smoothly cross over into the breathe- 
together solutions discussed in fToL fTT| . 



IV. SUMMARY 

In summary, we have derived an equation of motion for 
the total density fluctuations in dynamic two-component 
condensates. We have then used this to derive an ef- 
fective one-component CP equation for the smaller 
component, with effective potentials and interaction co- 
efficients which depend in a simple way on the relative 
scattering lengths in the system. We have studied and 
tested the predictions of this model in both effective re- 
pulsive and attractive interactions. Our model provides 
an intuitively simple way to predict the dynamics of two 
component BECs and allow us to make correspondences 
to results already obtained for one-component BECs. In 
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particular, our model provides new insight on the dy- 
namics of phase separation and the formation of vector 
solitons. We noted that magnetically trapped 87 Rb pro- 
vides a particularly interesting example in which the mo- 
tion is governed by a simple linear Schroedinger equation, 
allowing simple analytic predictions of motion in cases 
that would otherwise require solutions of the underlying 



nonlinear Schroedinger equations. Finally, we have ap- 
plied it to predict the existence of quasi-stationary two- 
component configurations. These solutions promise to 
be useful for applications involving information storage 
using stopped light as well as inducing spin squeezing in 
these systems. A combination of these two techniques 
may provide a new technique to generate squeezed light. 
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